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ABSTRACT 

The formation of planetesimals via gravitational instability of the dust layer in a protoplanetary 
disks demands that there be local patches where dust is concentrated by a factor of ~ a few xlO 3 
over the background value. Vortices in protoplanetary disks may concentrate dust to these values 
allowing them to be the nurseries of planetesimals. The concentration of dust in the cores of vortices 
increases the dust-gas ratio of the core compared to the background disk, creating a "heavy vortex." 
In this work, we show that these vortices are subject to an instability which we have called the heavy- 
core instability. Using Floquet theory, we show that this instability occurs in elliptical protoplanetary 
vortices when the gas-dust density of the core of the vortex is heavier than the ambient gas-dust density 
by a few tens of percent. The heavy-core instability grows very rapidly, with a growth timescale of a 
few vortex rotation periods. While the nonlinear evolution of this instability remains unknown, it will 
likely increase the velocity dispersion of the dust layer in the vortex because instability sets in well 
before sufficient dust can gather to form a protoplanetary seed. This instability may thus preclude 
vortices from being sites of planetesimal formation. 

Subject headings: accretion, accretion disks - hydrodynamics - instabilities planetary systems: for- 
mation - planetary systems: protoplanetary disks 



1. INTRODUCTION 

Current theories of planet formation postulate that 
dust grows from interstellar grain sizes (~ fim) to ~ km 
size planetesimals in the disks observed around young 
stars. This process must proceed in stages: below ~ 
cm scales, growth proceeds by "sticky" grain-grain col- 
lisions. Above ~ km scales, planetesimal growth pro- 
ceeds by gravitational accretion. However, around meter 
scales, collisions between grains are destructive and an- 
other growth process is needed. 

This process must be quick. The gaseous disk has a ra- 
dial pressure gradient which partially supports it against 
gravity, leaving its rotation rate sub-Keplerian. The 
dust, meanwhile, sees no pressure gradient and there- 
fore orbits at the Keplerian rate, leading to a headwind 
on the dust as it orbits in the gaseous disk. As a result, 
the dust rapidly spirals in to the star, at a rate 
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13/14 



yr 
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in the minimum mass solar nebula-far too short for 
the formation of plan ets (for a recent review, see 
IChiang fc Youdinll2009h . 

Gravitational instability in the dust layer is one 
mode by which pla netesim als can grow at the m scale 
(|Goldreich fe WardlflWl ISafronovl 11972) . However, for 
this to proceed the dust density must be significantly 
enhanced by a factor of ~ a few xlO 3 . The settling 
of dust to the midplane can enhance the dust den- 
sity enormously, but when the dust density, pd is sim- 
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ilar to the gas density, p g , K evin-Helmholtz instabili- 
ties limit further con centration ( Wcidcnschilli ng~fc Cuzzil 
119931: iChiand 120081 iBarrancd 12009ft . Hence, the dust 
density enhancement is limited to ~ 100 for solar metal- 
licity unless the me tallicity of the gas is enhanced (see 
lYoudin fc Shu! [20021 ) or the dust surface density is en- 
hanced. 

A natural way around these somewhat daunting 
timescale and surface density problems is to postulate 
the existence of regions in the disk where dust is sig- 
nificantly concentrated. This concentration may take 
place in persistent, large-scale vortices. On large enough 
scales (roughly > H, where H — c s /tt is the scale 

height of the disk, O = ^GM*/r 3 is the orbital fre- 
quency, M* is the mass of the star, r is the radius, 
and G is Newton's constant), the dynamics of a thin 
disk (H/r « 1) become quasi- two dimensional. Such 
flows host inverse cascade processes in which energy 
flows to large scales, making the appearance of coher- 
ent, long-lived vortices a distinct possibility. Such large 
scale vo rtices could be seeded by the baroclinic insta- 
bilities (iLpvelace et all 119991: iVarniere fc Taggen 120061: 
iLvra et alJl2009b iLesur fc Papaloizoull2009bj ). or bv the 
quasi-2D decay of initial turbulence generated from the 
initial ac cretion flow from th e pre-stellar envelope on to 
the disk (|Bracco et al.lll999D . 

Pre vious studies, both analytic (Barge fc Sommerial 
[l995t iTanga et all [l99l iChavanisI I2000D and com- 
putational dBracco et"aLl 119991 : iGodon fc Livid 120001 : 
ILvra et all 120091) . have demonstrated that anticyclonic 
vortices (those with vorticity antiparallel to the Keple- 
rian rotation) effectively trap dust. As these vortices col- 
lect dust, the dust-to-gas ratio in the vortices increases. 
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Hence, the gas in the vortices is denser than the gas in 
the surrounding disk. Can a density gradient between the 
vortex and surrounding gas or a density gradient within 
the vortex trigger new instabilities, destroying or mod- 
ifying the vortex in the process? As we will show in 
this paper, the answer is yes: density gradients in vor- 
tices are destabilizing for vortices sufficiently heavy cores 
(and also for all vortices with light cores). 

We note that this is not the only means by which 
dust can be concentrated into small regions of a proto- 
planetary disk. Youdi n fc Goodm an (20Q5|) have shown 
that the (inward) radial migration of dust in a gaseous 
disk is subject to a gas-dust streaming instability. The 
nonlinear evolution of this str eaming instability leads to 
large concentrations of d ust (|Youdin fe Johans"enl l2007t 
Uohansen fc Yo udin 2003), which might also be the sites 
of protoplanetary seed formation. As our primary in- 
terest is the stability of vortices, we do not study this 
mechanism, but mention it for completeness. 

This paper is organized as follows. We begin by dis- 
cussing equilibrium solutions for vorti c es in protoplane- 
tary disks in M focus ing on the iKidal (1981) ( 32~2|) and 
I Goodman et all (|1987L hereafter GNG) (t }2~T|) solutions. 
We calculate vortical stability in fj3j beginning with a 
description of Floquet theory ( H3.ip and then applying it 
to our equilibrium vortices in Wd.2\ We find two regions 
where vortices are unstable: vortices with light cores and 
vortices with sufficiently heavy cores. The growth rate 
of the instability is quite rapid for a sufficient density 
contrast. We discuss its application to protoplanetary 
vortices in 21 an d close with a summary of our results 
and discussion of outstanding issues in Sj5j 

2. VORTICES IN PROTOPLANETARY DISKS 

We model a local patch of the disk with a guid- 
ing center radius, ro, with angular velocity ft us- 
ing the incom pressible s hearin g sheet approximation 
(|Goldreich fc Lvnden-BelllH96l 
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(3) 
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where x = r—ro and y define the local coordinate system. 
We define the x and y velocities as u and v, respectively, 
and the gas pressure and density as P and p. Equations 
©, ©, (01), © are the continuity equation, x and y 
momentum equations, and condition of incompressibility, 
respectively. The following solution to equations ([!])-© 



u = 0, 
3 



(6) 

(7) 



defines the local backg round sheari ng flow. 

We now discuss the iKidal (|1981l) and GNG solutions 
to equations © - ([5]) . Both of these solutions have the 



v = —0JXX- 



(8) 

(9) 

Solutions that follow equation ((8} and ([9]) uniformly ro- 
tate on ellipses with ellipticity x >= 1 at an angular fre- 
quency, uj. For vortices in protoplanetary disks, u and 
ft have opposite signs: the vortices are anticyclonic. We 
note that while the Kida solution was originally derived 
for purely shearing flows in a fixed frame, i.e., for O = 
and without the tidal term, 3ft 2 a;, it has been applied to 
the study of vorti ces in protoplan etary disks as a means 
of collecting dust (]Chav anis 2000) and th e overa ll stabil- 
ity of vortices to 3-d effect (jLesur fc Papaloizou 2009a) . 
We solve for the equilibrium pressure distribution for 
the vortex solutions given by equations © and © by 
solving the steady state momentum equations © and 
Q. This gives 
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(10) 
(11) 



It is helpful to consider the pressure distribution in a 
coordinate system better suited to these vortices. As the 
steady state solution for both the Kida and GNG vortices 
have elliptical streamlines with ellipticity, x, we chose a 
coordinate system (6, <fi) of the form 

x = bcos<j), (12) 

y = &Xsin</>, (13) 

where b is the semi-minor axis that characterizes a par- 
ticular ellipse and <fr defines a position along that ellipse. 
Note that we have chosen our axes such that the y di- 
rection is the along major axis of the ellipse. In this 
coordinate system, 
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(14) 
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Using the above and equations (fTQ|) and (fTTj) . we find 
that 
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2.1. The GNG Solution 



GNG presented a solution for closed streamlines of the 
form given by equations ([8]) and ©. As this is simpler 
than the Kida solution (discussed below), we focus on 
this solution. Assuming a polytropic relation between 
the pressure, P and p, we find a relation between u>, ft, 
and X (GNG) 



ft 



X 



1 



(18) 



Applying the results of the above analysis of the pressure 
equilibrium (eq. [10], [11], US], and Q2]) to the GNG 
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vortex, we find: 

ldP 
p db 

1 dP 

bp d<j) 

The pressure distribution of the GNG vortex is very sim- 
ple compared to the Kida case. Its pressure gradient is 
zero along <p, and the pressure gradient is constant be- 
tween streamlines. Note, however, the pressure gradient 
in x-y coordinates is not constant moving along a stream- 
line due to their ellipticity. The pressure is negative out- 
ward (high pressure center) for x > 2 and inward (low 
pressure center) for \ < 2. Note that for x = 2, we find 
lo = — Q, which describes epicyclic motion and that the 
pressure gradient is zero, i.e., epicyclic motion demands 
no additional forces. 

2.2. The Kida Solution 



iKidal (fl98lh (see also IChavanisI l2000t 



iLesur &; Papaloizoul [2009a) presented an exact so- 
lution to the 2-D Euler equations (eqs.[5] - [5] with 
fi = 0) for a background shear and an elliptic patch with 
uniform vorticity, u). In the core, the vortex streamlines 
follow equations © and ©, while outside of the core, 
the streamline s asymptotic ally map onto the background 
shearing flow. IKidal (|1981l ) showed ( see IChavamsl 120001 
his Appendix A) that us, Xi an< 4 ^ f° r a time steady 
vortex is given by 



3n_ x( x -i) 

2w 1 + x 2 



(21) 



Fluid elements in these ellipses move at a constant an- 
gular velocity, 



ui 



ujx 



:m 



2(x-i)' 



(22) 



We should note that the IKidal (|1981f ) solution to the 
incompressible 2-D Euler equations enforces a nontriv- 
ial pr essure distribution (see also ILesur fc Papaloizoul 
l2009al) . Applying the same results of the above anal- 
ysis of the pressure equilibrium (eq. [10] , [TT] , [16] , and 
[17]) to the Kida vortex as we have done for the GNG 
vortex, we find: 



1 dP _ 3bCl 2 

P Ob 4(x-l) 

i op _ m 2 

bp d(j> 



of a terrestial vortex, where the base flow is axisymmet- 
ric and therefore allows us apply a cylindrical coordinate 
system. With this transformation the flow is simple and 
the perturbation equations separable. For spatially vary- 
ing flows such elliptical vortices, this is not possible and 
different techniques have to be brought to bear. 

One very pow e rful technique developed by 
iLifschitz fc Hameiril (]1991f ) combines Floquet anal- 
ysis with short wavelength WKB analysis. It is suitable 
for analyzing perturbations that grow both exponen- 
tially and algebraically in time or spatially varying flows. 
We provide a brief summary of this technique below 
and utilize it to analyze the stability of t he equilibrium 
vortex solutions of GNG and IKidal |l981l) . 

3.1. Floquet Theory 

H ere, we follow the logic o f Lifschitz fc Hame iri (1991) 
and iSipp fc Jacquinl (|2000D . We begin by perturbing 
inviscid incompressible Euler equations on the shearing 
sheet (eqs.[5] - [5]) to find 
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VP(27) 



where we have written the momentum equation in vec- 
torial format, and d/dt = d/dt + u • V. Now we will take 
perturbations of the form: 




exp 



i&(x,t) 




(x,t)+e 




(x,t) 
(28) 



where $ is a real phase function and e is a small param- 
eter. Inserting this ansatz into equation (|26p gives 



V$ u - 
V u 



0, 
= 0. 



(29) 
(30) 



to the lowest order and the next order in e, respectively. 
It is helpful to define the local wavevector, fc = V<&. 
Inserting the same ansatz into eauation (l25[) gives 






2 [(4X 2 - 8x + 7) cos 2 4> + 
(x 2 ~ 8x + 7) sin <\> cos 



3x 2 sin 2 , 



4x (X - #,) 



dp~ ~ „ 
-dt+ U - Vp 



0. 



(31) 
(32) 



The phase function $ is conserved along a streamline. 
Finato,) the perturbed momentum ()27[) gives 



4(X-1) 

The pressure distribution of the Kida vortex is rather 
complicated: the pressure gradient in the 4> direction is 
nontrivial, and for x > 4, the radial pressure gradient can 
vary from positive outward to negative outward where 
sin 0=1 (the long axis) . 

3. STABILITY OF PROTOPLANETARY VORTICES 

In general, spatially varying flows are not am enable 
to th e WKB-type analysis (see the discussion in iBavlvl 
1988). In the appendix, $A] we discuss a simple case 



du 
~dt 



u ■ Vk + 2f2 x u = — fc • P e 
P 



PP 
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Note that in deriving [33] we use the fact that to lowest 
order (1/e), the perturbed pressure gradient term in 1271 



gives 



ikP 

P 



0. 



(34) 



implying that P and thus VP are zero everywhere, as 
the local wavevector is in general not zero. 
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Fig. 1. — The basic vortex geometry. A sample streamline is in 
bold, with semi-minor axis b and semi-major axis c. The elliptic- 
ity X = c /^ = 4 in this case. <f>o shows the initial angle of the 
radius vector with the x axis while 4>k i s the initial angle of the 
perturbation wavevector k with respect to the radius vector. 

Now we project equation (j3"3")l onto the line perpendic- 
ular to k to eliminate P c using the operator / — kk/k 2 , 
where / is the identity. This yields 






■Zu- 



kk , 



' 20X "4(!f- / 



(35) 
where the velocity gradient tensor £ is defined as 

£ = Vu. (36) 

The other parts of the equation of motion are 

dp 



dt 



dx 

~dt~ 
dk 

~dt 



0, 



— = -£ T fc, 



(37) 
(38) 
(39) 



where we derived the last equation by applying V to 
equation (pTLj) . Equations ([33]) and (|3"T)) - ([59")) are the evo- 
lution equations for the position (x), perturbed density 
(p), perturbed velocity (u) and perturbation wavevector 
(fc) of the fluid perturbation. As x now has an explicit 
time dependence, this system of equations is amenable 
to Floquet analysis. 

3.2. Instability Analysis 

Taking the equilibrium solution (eq.[8] and [9]) posed 
above, we find 



£ = w 



X" 1 

-x o 



(40) 



Using this result to solve equations (|38|) and (|39j) . we 
find 



x = bcos(ujt + (fio), 
y = -Xbsm(ut + (j>o), 



(41) 

(42) 



k x = k cos(u}t + (f>k,o), 

k y = -x _1 fco sin(wi + <j> k ,o), 



(43) 
(44) 



where (f>o and <f>k,o are the initial phase of the coordinate 
and the wavevector and b is the semi-minor axis of the 
elliptical streamline (see figure [lj. We can also use the 
result k ■ u = (eq. [2H]) to write 

u x = a(b, i)sin(wi + 0fc,o), (45) 

u y =a(b, i)xcos(wi + </>k,o), (46) 

where a(b, t) determines the overall normalization of u. 
Applying the results of equations (|45|) - (|46|) to equations 
l(55|l and dHJ|, we find 

— = 2A~ 1 oj (x 2 - l) cos(ut + <j) kfi ) sm(ojt + (j>k,o)a 
+A _1 [(oj 2 + 3J7 2 ) cos(wi + <f>o) sin(u>t + <j>k.o) 
—X w " sin(w< + 0o ) cos(wt + 4>k,o) 



-2ttu>xam((f>k,o - 4>o)] b-, 
P 



dp 
~dt 



-asm((/) kfi - 4> ) 



dp 
db' 



(47) 
(48) 



where A = x 2 cos 2 (ujt + 4>k,o) + sin 2 (wi + <f>k,o)- Without 
loss of generality, we can choose 4>o = 0. Thus, equa- 
tions (|4"7|) and (|4^t depend only on the initial angle of 
the wavevector, <f>k,o- Equation (J4T|) explicitly depends 
on the semi-minor axis b. However, because b is inde- 
• Vijiendent of time and appears only in the p term in com- 
bination with p, we can eliminate it by absorbing it into 
the background density gradient, replacing dp/db with 
<91np/<91n& in equation ((48J) . 

The first term on the RHS of equation (|47|) results 
from the first two terms of equation (|35|) . The second 
term results from the density gradient within the vortex. 
For a zero density gradient, the integral over a period 
T = 2ircu~ 1 of equation (|47|) is zero as this first term on 
the RHS is an odd function. Hence, in the absence of a 
density gradient, no exponentially growing modes exists 
in two dimensions. In addition, for an initially aligned 
wavevector (i.e. <f>kfi — 0) and p(t = 0) = 0, there is 
also no growth, regardless of the presence of a density 
gradient. 

Equations (|48|) and (|47|) constitute a system of linear 
ODEs, which depends on initial conditions for a and p, 
the initial angle of the wavevector 4>k,o, the background 
density profile, dp/db, and the vortex ellipticity, x- How- 
ever, as we are interested in the asymptotic behavior of 
perturbations, i.e., growth or no growth, we are inter- 
ested in the parameter range of 4>k,o, Xi an d dp/db, which 
yield stability or instability. To this end, we define a 
growth rate 7 where pit = r) = e 7T /(r) and /(£) is a pe- 
riod function over the vortex rotation period, r = 2tt/uj. 
To find 7, we follow the general procedure of Floquet 
analysis on equations (|47|) and (|48)l . For linearly indepen- 
dent initial conditions, e.g., (a(t = 0) = l,p(t = 0) = 0) 
and (a(t = 0) = 0, pit = 0) = 1)0 we can integrate equa- 
tions (|47j) and (|48|) over one period 27r/w and solve for 

1 Note that we are ultimately interested in the growth factor 
over the period t, i.e., e 7T = p(r)/p(0) or = a(r)/a(0). Hence, it 
makes sense to rescale this linear problem in terms of the initial 
perturbation, i.e., a(t = 0) = 1 or p(t = 0) = 1. 
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Fig. 2. — Stability of protoplanetary vortices as a function of 
the density contrast. For any vortex that has a light core, i.e., 
dlnp/dlnb > 0, we find a purely growing instability, i.e. , the Vor- 
tical Rayleigh- Taylor Instability, which we discuss in i|3,3l For 
sufficiently heavy cores, there is an instability, which wc discuss in 
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Fig. 3. — Stability of protoplanetary vortices as a function of the 
X- For any vortex that has a light core, i.e., dlnp/dlnb > 0, we 
find a purely growing instability for any value of \ > 2. For suf- 
ficiently heavy cores, i.e., the <9mp/<91nb = —0.4 case, Instability 
also demands a sufficiently large \- 

the most unstable eigenvalue of the resulting matrix (e.g. 
iBender k Orszaglll978t iLesur k Papaloizoull2009a[ ). 

In Figure [5] we show the effect of different density con- 
trasts on the growth rate for fixed x = 10. Here we 
consider both light cores (d\np/d\nb > 0) and heavy 
cores (d\np/d\nb < 0). There are two regions of in- 
stability: one for light cores and and one for sufficiently 
heavy cores, which we discuss below. We show the ef- 
fect of the vortex ellipticity x on 7 m Figure [3] for the 
GNG (solid lines) and the Kida (dashed lines) for the 
light cores (light lines) and heavy cores (heavy lines). 
For the light core case, we find that as x increases, the 
growth rate in terms of u> also increases. However, there 
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Fig. 4. — Growth rates as a function of \ an d dlnp/dlnb for 
the GNG vortex. The heavy white contour denotes the stability 
boundary (y/ui = 0). 

is generally always growth (with the exception of the 
low x Kida case). In the heavy core case, there is no 
growth until x 1S sufficiently large and the core suffi- 
ciently dense compared to the ambient flow. We sum- 
marized the growth rate for (f>k.o — 7r/2 as a function 
of both x an d dlnp/dlnb in Figure ((H). In the case 
of the light core, d\np/d\nb > 0, we find instability 
for any density contrast. The heavy core case is more 
complicated. For small d\np/d\nb and small x> there 
is no instability. However, once d\np/d\nb and x are 
sufficiently large, growth sets in. Growth can occur for 
smaller d In p/d In b > if x is larger. In between the two 
regions of growth, i.e., light cores and sufficiently heavy 
cores, 7 < 0, indicating damping. The two heavy white 
lines marks the region of neutral stability. The physics 
of instability of these two light and heavy core cases are 
somewhat different, which we now discuss in £|3 .31 and 



3.3. Vortices with Light Cores 

The instability mechanism for light cores is analogous 
to the Rayleigh- Taylor instability. Recall that the pro- 
toplanetary vortices are high pressure regions. The pres- 
sure forces thus point outward and must be balance by a 
combination of centripetal and centrifugal forces which 
must point inward. If we " unroll" this vortex, we see that 
the combination of centripetal and centrifugal forces, 
which oppose the pressure force, is analogous to gravity 
in a pressure support atmosphere. Hence, a light core 
in this context is equivalent to making the material less 
dense where the pressure is largest, i.e., putting denser 
material on top of less dense material. This is subject 
to a Rayleigh- Taylor-like instability, which we call the 
Vortical Rayleigh- Taylor Instability (VRTI). In the ap- 
pendix, we present a simple example of this instability 
in a terrestial vortex to make more precise the analogy 
between the Rayleigh- Taylor Instability (RTI) and the 
VRTI. We note, however, that in the terrestial example 



() 
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presented in the appendix that the condition for instabil- 
ity is a heavy core as opposed to a light core. This results 
from the fact that terrestial vortices are low pressure re- 
gions, whereas protoplanetary vortices are high pressure 
regions. 

We now demonstrate the behavior of perturbations for 
light cores by plugging to for GNG (or Kida) vortices 
(eq. [18]) into equation (|47|) . we integrate the evolu- 
tion of equations (|47|) and (|48|) over several periods. We 
note that equation (|47[) admits a purely analytic solution 
when p = 0. Integrating both sides with this in mind, 
we find 



a(t) = a 



X 2 + (1 - X 2 ) sin 2 (fo, ) 
X a + (l-X 2 )sin 2 (wt + &,o)' 



(49) 



where ao is fixed by initial conditions. The analytic solu- 
tion (|49|) represents the evolution of a perturbation that 
is purely advected along in the flow. More complex cases 
are solved numerically. 

Figure [5] shows the behavior of a(t) and p(t) for the 
initial conditions a(t = 0) = 1 and p(t = 0) = and a 
background density profile of d In pj d In b — 0.01, i.e., a 
light core. The background density profile corresponds 
to a one percent decrease in the density across the vor- 
tex. The velocity amplitude a(t) for the 4>k,o = case 
shows oscillatory behavior between 1 and 100, but no 
long term growth occurs. Correspondingly, the density 
perturbation remains zero for all time and is not shown 
in the right panel. Indeed, the solution's behavior pre- 
cisely follows the analytic result of equation (j49|). veri- 
fying the accuracy of our numerical integration. On the 
other hand, both 4>k.o = 7r/4 and <ftk,o = tt/2 show signif- 
icant growth, and the asymptotic growth rate, i.e., the 
slope of the trend, is maximized for n/2. This is un- 
surprising given the discussion in the appendix and the 
analogy with the RTI, where we found that growth rates 
are maximized when the wavevector is parallel to vortex 
streamlines. 

Similarly, we plug in uj for Kida vortices (eq. |22| ) 
into equation (|47[) and integrate its evolution over sev- 
eral periods. Before comparing the behavior of the Kida 
vortices with that of the GNG vortices, we first note that 
from Figure [5l the velocity and density shows both short 
timescale periodic fluctuations, a result of the spatially 
inhomogeneous flow field, i.e., vortex, and long term be- 
havior. As we are only interested in long term behavior, 
we sample p where tot is an integer multiple of 27r for both 
vortices and show the results in Figure |6l The compar- 
ison clearly shows that the background state makes no 
difference in the qualitative behavior of the instability 
and little difference in the quantitative behavior of the 
growth rate. 

3.4. Vortices with Heavy Cores 

We now discuss the heavy core case. To help eluci- 
date the physics, we first consider the effective gravity 
that counteract the pressure forces in a GNG vortex, 
i.e., p~ 1 dP/dx — —g x , p~ 1 dP/dy = — g y . From equa- 
tion (|19[) , we know that the pressure drop between vortex 
streamlines is constant. Hence, the effective gravity is 



^dPdb 
~ P ~dbdx~ 



, dP db 



(51) 



where b = y ' x 1 + y 2 j\ 2 . Hence computing the magni- 
tude of the effective gravity is, thus, 



■9v =9o\ /cos 2 . 



where 



go = bn 2 



3X 2 



x 2 -i 



-2« 



3x 2 



x 2 -i 



= const. 



(52) 



(53) 



(50) 



We now plot the behavior of the effective gravity 
(§(t) = g(t)/go) and the growing perturbations (a(t) and 
p(t)) in Figure [7] for 4>k,Q — 7r /2- Note that perturbed 
density {pit)) and velocity [ait)) changes only at inter- 
vals where g(t) is large. This is unsurprising for large 
X, the effective gravity will vary between go/x an d So- 
These peaks in the effective gravity is reached for y = 0, 
i.e., on the minor axis of the elliptical vortex. This is rea- 
sonable from a physical perspective as it is here that the 
distance between vortex streamlines is minimal while the 
pressure change between vortex streamlines remain un- 
changed (for GNG vortices). Thus, the force is maximal 
there. This effective gravity, which is time dependent (as 
fluid is advected along a vortex), is akin to a kick. The 
relative strength of these kicks depend on \ an d the pe- 
riod of these kicks is exactly half of the rotation period of 
the vortex as a fluid element moving along a streamline 
crosses the minor axis twice a rotation period. As the 
modes we are following is akin to radial gravity modes, 
the density contrast determines their period. Hence, the 
minimum density contrast required for instability has an 
obvious intepretation: the radial gravity mode must also 
have a period that is comparable to the rotation period 
of the vortex to it to couple successfully to the kick and 
develop overstable oscillations as the plot of p in Figure[7] 
shows. In addition, the minimum \ needed for instability 
as demonstrated in Figure [3] and |4] results from requiring 
each kick to be sufficiently strong. 

Similarly, we can make the same comparison between 
the GNG and Kida vortex in the heavy core case in Fig- 
ure [5] as we have done in the light core case in Figure 
[6] Again, we plug in to for Kida vortices (eq. [22]) into 
equation (|47|) and integrate its evolution over several pe- 
riods but this time for a background density gradient of 
d\np/d\nb — —0.3 and \ — 12. Again we discard short 
timescale periodic fluctuations and sample p where tot is 
an integer multiple of 27r for both vortices. The com- 
parison shows that the growth rate of the HCI is more 
dependent on the vortex solution (Kida vs. GNG) than 
the VRTI. Whereas Figure [6] shows that the amplitude of 
the density perturbation for the Kida and GNG vortices 
track each other fairly closely, these amplitude diverge 
much more strongly in the HCI as shown in Figure [5] 

4. DISCUSSION 

The analysis of the preceding section demonstrates 
that vortices with light cores or sufficiently heavy cores 
are unstable to the VRTI and HCI respectively. Figures 
[2] and [4] shows that growth occurs on a few vortex ro- 
tation periods. Moreover, these instability appears to 
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Fig. 5. — Velocity (left) and density (right) perturbations as a function of time for three initial wavevector orientation angles (f>f. q for 
light cores (91np/91n6 = 0.01). Note that <j>k = h as been dropped from the right plot because it has zero growth. This case is show on 
the left plot for illustrativ purposes. There is no growth for the 0^ q = case, and the growth rate is maximum for for <f>k = Jr/2, 
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Fig. 6. — A comparison of the GNG and Kida background states 
for the <f>k q = t/2, \ = 5 case for the light core case (d In p/d In b = 
0.01). The figure shows density perturbation as a function of time 
for both background states. 

be robust and its detailed physics are independent of 
the vortex model used (either GNG or Kida). Having 
demonstrated the basic physics of these instabilities, we 
turn now to its application to planetesimal formation. 
We will first review some of the physics of dust trapping 
in vortices. 

If planetesimals form by gravitatio nal collap s e and 
fragmentation of a du st sublayer (jSafronovl Il972t 
IGoldreich fc Ward! I1973J ). then this layer m ust have a 
Toomre, Q < 1 (but also see iWarcll 12000ft . which im- 
plies that the velocity dispersion of the sublayer must be 
below: 

.^^^cms- 1 (54) 

for the minimal mass solar nebula (MMSN), where Sd 




2.0 2.5 

t -u/2-k 

Fig. 7. — The effective gravitational acceleration g(t) = g(t)/go 
normalized to the maximum gravitational acceleration, perturbed 
density (in arbitrary units), and a(t) (again in arbitary units). Note 
that the peaks of g(t) correspond to the peaks in a(t) and changes in 
p(t), which correspond to periodic kicks at roughly half the rotation 
period of the vortex. 

is the surface density of the dust layer. For a lam- 
inar disk, this criterion is amply fulfilled if the dust 
is allowed to settle to the midplane. However, as the 
dust collects near the midplane, it is subject the in- 
duced Kelvi n-Helmholtz instabilities with t he overlying, 
gas layers (iWeiden schilling fc Cuzzl 119931 : iCuzzi et"aLl 
HL993I : lChiangl r2008: Ba rrancdl2009T r Therefore, the dis- 
persion of the dust layer is closer to a few ms _1 . Hence 
the surface density must be enhanced by a factor of 
~ 20 — 100 (effectively the Q of the g aseous disk) so t hat 
gravitational instability can operate ()Chavani s 2000) . 

More careful considerations suggest this enhancement 
of ~ 20 — 100 may be a severe overestimate and that 
only enhancement of order a few is needed in the 
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Fig. 8. — Same as Figure [6] except for the heavy core case. Here 
X = 12 and 9 In p/d In 6 = -0.3. 



high metallicity disks which preferen tially form planets 
(jYoudin fc Shull2002t [Johansen et al.ll2009Tl . In any case, 
vortices are one avenue by suc h a dust surface density en - 
hancement can be achieved as lBarge fc Sommerial (1995) 
first pointed out. The timescale for dust to connect and 
concentrate in vortices, i.,e., the capture timescale, i C apt, 
can be fairly ra pid, i.e., £ caDt ~ t Awni when t s tOD ~ t Awn as 
pointe d out by iBarge fc Sommerial (|1995[ ); iTanga et al.l 
i|1996tl : iChavanisI ([20000 . Over "the lifetime of a vortex, 
iiife, the amo unt of dust the c an be gathered by a vortex 
is very large. [Chavanisj ([20000 argues that this mass is 



Aid ~ ^^lifcSd-R / (istop), 



(55) 



where / describes the efficiency of capturing dust and 
is ~ 1 when i s t op ~ idyn and R is the size scale of the 
vortex. If the inward concentration of dust is balanced 
by the outward diffusion of this dust concentration due 
to turbulence, these dust particles would be confined to 
a region on a scale of 



r d ~y/Dt 



capti 



(56) 



where D ~ aR 2 fl is the turbulent diffusivity. For 
a = 0.01, rd ~ O.li?, i.e., in the central core of the vor- 
tex. Hence the surface density is enhanced by two orders 
of magnitude. Since the GI hypothesis for planetesimal 
formation only demands more modest increases, vortices 
should be ideal sites of planetesimal formation. 

However, such a increase in dust surface density is not 
without its costs. As we have shown a sufficient increase 
in the effective mean molecular weight of the gas in the 
cores of vortices is destabilizing. The condition for the 
HCI demands a mean molecular weight increase of or- 
der 20% or an increase of the dust surface density by a 
factor of 20% if the dust and gas densities are similar 
in the midplane. This increase in dust surface density is 
much smaller than what is required for the GI hypothesis, 
which demands a factor of a f ew increase (jYoudin fc Shul 
120021: iChiang fc Youdinll2009Q if the vertical structure of 
dust is taken into account to a factor of 20 — 100 when 
vertical structure is ignored ()Chavanisll2000D . Thus, the 



HCI will be triggered before gravitational instability sets 
in according to the present linear calculation. 

There are many issues involving the stability of vortices 
with a heavy core than cannot be resolved by the present 
linear calculation, which we now briefly discuss. The first 
issue is the non-linear state of the instability, which is not 
known at present. We expect the HCI to grow until its 
saturates, which may 1. destroy the vortex, 2. increase 
the velocity dispersion of the dust layer, or 3. limit the 
enhancement in the dust surface density to a few tens of 
percent, i.e., marginal stability. For any of these options, 
GI is curtailed in cores of vortices. 

Another issue is the equilibrium distribution of dust 
along a streamline. We have assumed that the dust is 
uniformly distributed along a s treamline. For light par- 
ticles, this is li kely the case. ITanga et al.l (1996) and 
IChavanisI ([20000 studied the process of dust trapping in 
vortices and found that the zeroth order motion is that 
light particles of dust travels along the elliptical stream- 
lines with a slow "radial" drift due to drag forces. How- 
ever, heavy particles move along epicycles, i.e. , ellipses 
with aspects ratio 2. In addition, Youdin (2008) showed 
that the stationary point for dust in a sub-Keplerian gas 
is not the center of the vortex, but rather a point that 
is forward in azimuth. This is unsurprising as the dust, 
in maintaining a sub-Keplerian rotation rate, demands 
an additional radial force away from the central star to 
counteract gravity, a force that is supplied by gas push- 
ing on the dust if the dust is ahead of the vortex center in 
azimuth. These elements suggest that the distribution of 
dust along a streamline may not be uniform and so may 
affect the stability properties of vortices in a non-trivial 
way. 

A third issue is the nature of gas-dust coupling. We 
have assumed the gas and dust are well coupled on a 
dynamical time. However, this may not be the case. 
For instance, the fastest settling dust is that which is 
marginally coupled to the gas, i.e., H,t s = 1, w here t s 
is the dust stopping time, (Johans en et al.l [20 041. The 
dust that is trapped in vortices may be preferentially of 
a certain size, i.e., marginally coupled to the gas. Hence 
the instability growth time, dynamical time, and dust- 
gas coupling time, in the dusty protoplanetary disk can 
be all of the same order. 

Additional instabilities that arise directly from 
this gas-dust coupling may also be important. 
lYoudin fc Goodman (2005) showed that the imperfect 
coupling between gas and dust and their backreaction 
on each other leads to a secular streaming instability. 
This instability leads to protoplane tary disk turbulence 
and te nds to concentrate dust lYoudin fc Johansenl 
(|2007fk Uohansen fc Youdin! ([2007Q . These streaming 
instabilities or analogues may also be important in the 
stability protoplanetary vortices and would be profitable 
to explore. A proper accounting of gas-dust coupling in 
a vortex and its effect on the VRTI is a topic of future 
work. 

Another issue is that the vertical structure of the gas 
and the dust is very different in protoplanetary disks. 
Dust tends to settle toward the midplane (though the 
presence of vortices and/or turbulence may counter this 
tendency). This settled dust may drive vertical turbu- 
lence if it is sufficiently concentrated (see for instance 
iChiangj 120081; iBarrancd [20091) . The effect of this differ- 
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ence in the vertical structure of gas and dust on vortices 
has not yet been addressed and is likely important for 
both their equilibrium and stability. However, we may 
expect the HCI to be important regardless because is a 
2-D effect in a thin dust layer within a thicker gas vortex. 
Finally, the alert reader (and referee) will note that we 
have not discussed the case of 1 < x < 2 vortices, i.e., 
low pressure vortices. In principle, such vortices might 
exist in protoplanetary discs, but the prevailing theoret- 
ical bias is that vortices in protoplanetary discs are high 
pressure regions. We have gone along with this bias in 
this work and have ignored these low pressure vortices. 
However, we note that heavy core low pressure vortices 
are violently unstable (equivalent to the light core case 
discussed above) because of the reversal in the direction 
of the effective gravity. In addition, it is unclear if these 
vortices would concentrate dust. The equivalent calcula- 
tion of lChavanisI (J200Q[ ) for low pressure vortices has not 
been performed. While a detailed study of low pressure 
protoplanetary vortices might be interesting, the impact 
of such a study is unclear. 

5. CONCLUSIONS AND OPEN ISSUES 

We have demonstrated two instabilities in protoplane- 
tary vortices, resulting from light cores (VRTI) and suf- 
ficiently heavy cores (HCI). The physics of the VRTI is 
analogous to the Ray leigh- Taylor instability, with gravity 
replaced by centrifugal and centripetal forces in a rotat- 
ing fluid. The HCI appears to be a parametric instability. 
We have shown that these instabilities are robust for all 
vortices possessing a light or sufficiently heavy core. For 
protoplanetary vortices, the instability of interest is the 
HCI as dust would concentrate in their centers, leading 
to heavy cores. While the nonlinear state of the these 
remains unexplored, we expect that this instability pre- 
vents vortices from acting as protoplanetary nurseries. 

Both the VRTI and HCI are novel among elliptical 
vortex instabilities as they are 2-D - only motions in 
the x-y plane are required. Previous work on the stabil- 
ity of vortices have focused on the importance of the 



instabilities that involve 3 -D motions (|Lithwickl 120091 : 
iLesur fc Papaloizoull2009a|) . Indeed, 3-D effects do lead 
to additional instabilities t hat destroy vortices even be- 
fore they can collect dust. ILesur fc Papaloizoul ((2009a) 
argued that the 3-D elliptical instability can destroy Kida 
vortices for x < 4 and x > 6. Lithwick (2009) argues that 
nonlinear coupling and transient amplification between a 
vortex and its children that involve a vertical component 
leads to destruction of the vortex. He argues that the 
stability requirement for any vortex is then that the base 
of the vortex be larger (by at least a factor of 2) than 
its height. Our results s uggest tha t if gas vortices having 
X > 6 (as proposed by iLithwickl (120091 )). the HCI will 
affect these vortices as they gather dust. 

The analysis that we have attempted here is linear and 
so it is highly dependent on the background equilibrium 
state. The two equilibria (Kida and GNG) that we have 
analyzed in this paper were chosen due to their simple 
analytic structure. Although we have shown that the 
HCI is very similar in both these cases, 3-D simulations 
have clearly demonstrated that vortices in protoplane - 
tary disks are not so sim ple Bar ranco fc Marcus! (2005); 
iShen et all (12006ft : iLiThwickl (l2009f l~ 

Finally, our work leaves open a number of issues in- 
cluding the nonlinear state of the HCI, gas-dust coupling 
physics, and equilibrium structure and vertical structure 
of vortices. We are currently pursuing numerical work 
exploring the effects of heavy vortices in protoplanetary 
disk. 
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APPENDIX 
A SIMPLE EXAMPLE OF VORTICAL RAYLEIGH- TAYLOR INSTABILITY 

In this appendix, we discuss a simple example of the Vortical Ray leig h- Taylor instabilit y (VRTI) in a terrestrial 
vortex to illustrate its basic physics. Our simple treatment is derived from lSipp et al.l (|2005l ). and a detailed overview 
of the state of terrestrial heavy vortex instability theory is found therein. 

For a 2-d circular vortex in equilibrium, pressure forces (inward) are counterbalanced by centrifugal forces (outward). 
Incompressible motions of this 2-d vortex are described by the continuity equation, 



dp dp v dp 

— + u— + -— =0, 

dt or r do 



(AI) 



where p is the density, u = r and v — r<j), the momentum equations, 



dt dr r deft r p dr 



dv dv v dv 
dt dr r dd> 



uv I dP 

r pr d(f> ' 

where P is the pressure, and the incompressibility condition, 

i dru 1 dv 



dr 



r dm 



(A2) 
(A3) 

(A4) 
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Fig. 9. — A cartoon showing the salient feature of the Rayleigh-Taylor type instability, which is that the density gradient (represented 
by the shading: dark tones correspond to high density, light tones to low density) is in the opposite direction from the restoring force. In 
the case of the standard Rayleigh-Taylor instability, the force is gravity; in the vortical Rayleigh-Taylor case it is the centrifugal force. 

In equilibrium, the fluid motions of the vortex are circular and constant, i.e., u — and v is constant. Hence we find 
that 



dr 



v 
r 



We now perturb equations (IA1|) - (IA4|) and assume perturbations of the form exp(— iujt + im(j) 
continuity and momentum equations read 



iaSp - 
-iaSu 



dp 
dr 



= 



2v 



-Sv ■■ 



-ik 



SP 



ldPSp 
p dr p 



. dv v \ .771 SP 

-laSv + \ — — I — ou= —i , 

dr r J r p 



(A5) 
ikr). The perturbed 

(A6) 

(A7) 

(A8) 



where a — (u> — mv/r). The incompressibility condition feq. |A4p becomes 

■i r * m c 

ikou H ov = U, 

r 

where we have assum ed fcr » 1. Equation (|A9j) gives Sv in terms of Su, which we apply to equations (|A7|) and (|A8 
Using equation (|A6[) for 5p in terms of Su, we find the dispersion relation: 



(A9) 



m 



iak- 



dr 



m v <91n p 
r 2 r dr 



0. 



(A10) 



For an uniformly rotating vortex, v oc r, the second term in (jAlOj) vanishes and the solution to the dispersion relation 



is 



a 2 = 



m 2 /r 2 v 2 dlnp 



• ■ , , -, • ( Al1 ) 

m z [r z + k z r or 

which is < (unstable) if dlnp/dr < 0, that is, if the core of the vortex is heavy0 This instability is analogous to 
the Rayleigh-Taylor instability, whose dispersion relation is w oc k\jk 2 : , but where gravity is replaced by a centrifugal 
force. In the case of the Rayleigh-Taylor instability, the equilibrium is set by pressure forces balancing gravity. Heavy 
fluid that sits on top of light fluid which fulfills the conditions of equilibrium, but is unstable to interpenetration across 
the interface. By analogy, in the VRTI case, the equilibrium vortex is set by pressure forces balancing centrifugal 
forces. The presence of a heavy core leads to non-axisymmetric instabilities (where the origin is set by the center 
of the vortex), where again the heavy fluid elements in the core interpenetrate light fluid on the exterior. Figure [S] 
illustrates this analogy. 

The analogy is made clearer if we identify k^ = m/r, which is perpendicular to the centrifugal force (i.e., the effective 
gravity), which is in the f direction. Thus, we can make the identification k\j(k\ + k 2 ) <H> fc^/fc 2 between VRTI and 
RTI, respectively. The wavevector that grows the fastest in both instabilities is the wavevector that is perpendicular 
to the vertical gravity and the radially outward centrifugal force, respectively. 
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